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. . ABSTRACT 

O 

O 

Tj- ' The interaction of a stellar or disk wind with a collapsing environment holds 

promise for explaining a variety of outflow phenomena observed around young 

stars. In this paper we present the first simulations of these interactions. The 

t;;-!- ' focus here is on exploring how ram pressure balance between wind and ambient 

^ ■ gas and post-shock cooling affects the shape of the resulting outflows. In our 

O ' models we explore the role of ram pressure and cooling by holding the wind 

-J. ■ speed constant and adjusting the ratio of the inflow mass flux to the wind mass 

0\ • flux {Ma/ Myj) Assuming non-spherical cloud collapse, we find that relatively 

^ ■ strong winds can carve out wide, conical outflow cavities and that relatively weak 

I ' winds can be strongly coUimated into jet-like structures. If the winds become 

;_i ■ weak enough, they can be cut off entirely by the infalling environment. We 

^ ■ identify discrepancies between results from standard snowplow models and those 

cd ■ 

presented here that have important implications for molecular outflows. We also 



present mass vs. velocity curves for comparison with observations. 

Subject headings: ISM: infall and outflows - ISM:molecules - stars: young - stars: 
winds, wide angle 



1. Introduction 

Molecular outflows are commonly observed in association with young stars. Although 
the precise mechanism generating the outflows is poorly understood, it is generally believed 
that the molecular mass is driven and excited into emission by a wind emanating from the 
inner circumstellar disk close to the star, and driven by accretion energy (Calvet & GuUbring 
1998). Recent observations have explored how these outflows form close to the driving source 
(Chandler et al. 1996) and how they might evolve in time (Velusamy & Langer 1998). These 
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and other investigations indicate that molecular outflows are not clearly explained in terms 
of any single current model. 

Molecular outflow models are principally divided into two groups; those in which out- 
flows are driven by highly-collimated winds or "jets", and those in which outflows are driven 
by wide-angle winds. Masson & Chernin (1992, 1993) and Chernin & Masson (1995) com- 
pared jet driven and wind driven models of molecular outflows. They conclude that jet-driven 
models are a better flt than wide-angle wind models because observations show very little 
outflow mass moving at the highest velocities. A jet with moderate post-shock cooling would 
transfer its momentum to the ambient molecular material largely through the small cross- 
section of its head, while a wide wind has a much larger driving area. On the other hand, 
jet-driven models have so far had some difficulty producing the cross-sectional widths of many 
observed outflows. To reproduce these morphologies, investigators have invoked mechanisms 
such as turbulent entrainment (Stabler 1994; Raga & Cabrit 1993), and precession (Cliffe 
et al. 1995; Suttner et al. 1997). These strategies imply fllled outflow cavities and shock 
structures that aren't generally observed (Chandler et al. 1996). Wide angle wind models 
can produce wind-blown bubbles with wide bow shocks and empty cavities. It is, however, 
difficult for the simplest versions of these models to reproduce the observed momentum dis- 
tribution. Li & Shu (1996) have suggested that a wide angle wind whose properties vary as 
a function of polar angle may help solve the momentum problem. Such a wind requires a 
structure that has higher density along an axis: more like a jet. Thus at present, models 
and observations leave it unclear as to what is really the driving mechanism of molecular 
outflows. With these issues in mind, the present paper examines a wide angle wind model 
through hydro dynamic simulations. 

The early model used by Chernin and Masson to make their case against wide angle 
winds is that presented by Shu et al. (1991). It applies a "snowplow" scenario where a central 
wind sweeps up ambient material and compresses it into a thin shell. The calculation is 
performed under the critical assumptions of isothermal shock dynamics and strong mixing 
between the post-shock wind and ambient gas. The model successfully produces collimated 
bipolar wind-blown bubbles when relatively denser "equatorial" regions prevent the shell 
from expanding at low latitudes as quickly as it does along the pole. 

New features have been added to the basic wide angle-wind blown bubble model. Close 
to the young star, effects such as gravity from the central star and the density distribution of 
the inner envelope become crucial in influencing the shape and evolution of the outflow. This 
has been investigated in detail by Wilkin & Stabler (1998), who model a quasi-static balance 
between a steady central wind, gravity, and a time-dependent angularly varying infall. In 
essence the bubble is considered as a series of accretion shocks. Using the momentum balance 
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across the fixed shell face in both normal and transverse directions, the authors solve for 
the density and flow pattern within the shell. Because of the quasi-static assumption, the 
time-dependence of the bubble shape and size is directly linked to the way in which the 
infall changes in time. The stability of such a situation remains to be investigated. While 
they found that signiflcant collimation can occur as the environmental density distribution 
becomes more flattened, the timescale for this is ~ 10^ yr. They point out that this is much 
longer than what is observed. This result reveals the importance of performing dynamic, 
time-dependent calculations. 

While considerable progress has been made in the wide angle wind scenario for molecular 
outflows there remain many aspects which have not yet been explored. In particular the 
full time-dependent mult i- dimensional nature of the flows has not been examined. Solving 
cylindrically symmetric non-linear, time- dependent fluid equations by numerical simulation 
we obtain signiflcantly different behavior from that obtained in previous wide angle wind 
blown bubble models. The shape of the envelope, namely the gradient of density with polar 
angle, not only affects the shape of the outer shock and resulting molecular shell, but also 
the shape of the inner (wind) shock. This can have important dynamical consequences as 
has been shown in Frank & Mellema (1996) and Mellema & Frank (1997). 

It seems difficult to avoid the the conclusion that some intrinsic collimation of the initial 
mass ejection is required, given the evidence for highly-collimated jets over long scales (Bally 
et al. 1997; Bachiller 1996). This is particularly true for emission-line jets of optically- visible 
T Tauri stars (e.g., Stapelfeldt et al. 1997; see also the review by Reipurth & Heathcote 
1997), where there is little evidence for the presence of enough external (dusty) medium for 
initial hydrodynamic collimation of ionized jets to be effective. On the other hand, it should 
also be recognized that even models of mass ejection with collimated flows often have some 
wide-angle component. This is true not only of the magnetic "X-wind" model of Shu and 
collaborators (Shu et al. 1994; Najita & Shu 1994; Shang et al. 1998), in which a dense axial 
flow is surrounded by a lower-density dispersing flow; it is also true of the original self-similar 
MHD disk wind (Blandford & Payne 1982; Pudritz & Norman 1983; Kdnigl 1989). For this 
reason it is worth exploring the interaction of wide-angle flows with the ambient medium. 

In this paper we consider the time- dependent interaction of an infalling envelope with 
an expanding wind. Our focus is spherically-symmetric winds as an example of a maximally 
wide angle wind, though we do calculate two cases using aspherical winds for comparison. 
Collimation of the resulting flow will be enhanced if the infalling envelope is not spherically- 
symmetric. We take one of the simplest models for such non-spherical collapse, the model 
of Hartmann et al. (1996), from an initial self-gravitating sheet with no magnetic fleld. 

In § 2 we present the physics and the computational scheme used for the simulations. 
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The simulation results are examined in § 3. In § 4 we compare our simulations to the 
snowplow model of outflows, discuss the possible role and importance of turbulence in these 
models, and make comparisons of the simulations to observations. We conclude in § 5. 



2. Computation 

2.1. Methods 

The simulation code solves the Euler ideal fluid equations with source (sink) terms for 
central point-source gravity and radiative cooling. The equations, written in conservative 
form, are: 



§^ + V-pu = 0, (1) 

^ + V ■ puu = pV<l>, (2) 

^ + V ■ u(e + p) = -A(p, T) - pu ■ V$. (3) 

where p is the fluid mass density at a point, and u is the velocity. The pressure p is related to 
the energy density e and its kinetic (efc„) and thermal components (eth) through the relations 

e = Cfen + eth, (4) 

Cfcn = -PM^ (5) 

The adiabatic index 7 = 5/3 is that for a monatomic gas. The temperature is determined 
through the ideal gas equation of state, with the particle mass set to that of atomic hydrogen: 

pkT 
p= . 7 

The cooling source term A(p, T) includes a Dalgarno-McCray 1972 cooling curve which is 
applied above 6000 K, and a Lepp and Shull 1983 cooling curve applied below 6000 K to 
simulate the effects of both high and low temperature cooling. 

Gravity in these models is purely due to a source at the origin. Self-gravity of the fluid is 
not included. The gravity source term is written in terms of the gradient of the potential 
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V$. This force is specified via the constant central source mass M* as 

V$ = G^f, (8) 

The radius r is the vector from the origin to the fluid point. 

The equations are solved on an Eulerian grid using an operator split numerical method, 
where several operators are applied per timestep, each one simulating a different aspect 
of the physics. We also employ fluid tracking to determine what part of the initial grid 
a fluid parcel came from. The grid covers a quarter meridional plane of a cylinder and 
has assumed axial symmetry plus reflective symmetry across an equatorial plane. The basic 
hydrodynamics equations without source terms are applied via the numerical Total Variation 
Diminishing (TVD) method of Harten 1983 as implemented in dimensionally split form by 
Ryu et al. (1995). The source term for radiative cooling is applied explicitly to first order 
via an exponential as described in Mellema & Frank (1997) where: 



er=elkexp ^^^^At . (9) 



'th 



The superscripts represent a value at the nth and (n + l)th timestep, and At represents 
the amount of time between timestep n and n + 1. Gravity is applied through a first order 
explicit Euler operator which updates the velocities and the kinetic energy density. The fluid 
is evolved by rotational velocity through another first order explicit Euler operator. 

The timesteps are adjusted to satisfy the Courant condition on the sourceless hydro- 
dynamics, and simultaneously to be less than a reasonable multiple of the cooling time: 
3A(p, T)/et/i. Because of the high densities in the ambient medium [n > 10^ cm^^) the 
cooling time step can become relatively small and the number of times steps required to 
complete the simulation relatively large. This produces additional diffusion in the solutions 
causing our flows to appear rather smooth. By artificially reducing the cooling in some test 
simulations, we were able to obtain sharper features. The gross morphology and evolution 
of the flow however remained qualitatively the same as those with the full cooling shown in 
§ 3. 

We applied several tests to the code. We compared momentum conserving wind blown 
bubbles with the solutions obtained by Koo & McKee (1992) to test the basic hydrodynamics 
with strong cooling. We performed radiative shock tube tests along the axis using a simplified 
power law cooling function and compared these against semi-analytic shock solutions to test 
the cooling operator. To test gravity we reproduced the Bondi accretion solution in the code. 
To test gravity and rotational velocity, we maintained a Keplerian ring in an orbit around a 
central mass. In general, the code was able to recover the analytical solutions to better than 
10%. 
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2.2. Model 



The scenario we applied involved driving a spherically symmetric wind into an infalling 
non-spherical envelope. While a spherical wind may not be operating in many young stellar 
object systems, we seek in this investigation to explore the shaping of the outflow through 
wind/environment interactions and make contact with previous analytic studies. We hope to 
relax the assumption of a spherical wind in a later paper. The parameters used are given in 
Table 1 and are based on typical values for low mass young stellar objects. A radially directed 
steady wind of velocity 200km/s and mass flux 10^'^ MQ/yr was imposed at each timestep 
on the cells of the grid within an origin-centered sphere of radius equal to 10% the radial 
size of the simulation. We chose an environment given by Hartmann et al. (1996) derived 
from their simulations of a collapsing, rotating, axially symmetric sheet. In this model, the 
collapse of the initially-flattened protostellar cloud gives rise to a highly anisotropic density 
distribution in the infalling envelope, which in turn helps strongly collimate the initially 
spherical flow. As Hartmann et al. point out, whether or not this model is correct in 
detail, the general property of non-spherical clouds to flatten or become more anisotropic as 
gravitational collapse proceeds lends credence to the idea of aspherical infall. The anisotropic 
stresses of magnetic flelds can also produce non-spherical protostellar clouds, and Li & Shu 
(1996) argue that such initial conflgurations can also help make outflows more collimated. ^ 

The density is given by the equations (8), (9) and (10) in Hartmann et al. (1996) which 
modifies calculations by Cassen & Moosman (1981) and Terebey et al. (1984) (referred to 
hereafter as TSC). The velocity distribution is taken from Ulrich (1976). The equations 
are rewritten in a slightly different but equivalent form in the appendix. These equations 
produce a flattened infalling toroidal density distribution with an equator to pole density 
contrast Pe/Pp > 1000. This is larger than what is produced by the Cassen & Moosman and 
TSC models. The TSC model was used by Wilkin & Stabler (1998). 

The environment near the center of the sheet models used in this paper is a torus of a 
high equatorial density (6 x 10^ iriH/cm?) with a half-maximum opening angle of almost 180 
degrees. In the hydrodynamic coUimation simulations of Frank & Mellema and Mellema & 
Frank a "fat" torus was used with opening angles of ~ 90°, similar to those obtained by Li 
& Shu (1996) for magnetized collapse. These previous simulations produced outflows with 
strong collimation. It is noteworthy that, as we shall see in § 3, it is possible to produce 



^However, we note that the static configurations suggested by Li & Shu may not be particularly relevant 
to outflow sources, because some collapse must have already occurred to produce the central mass responsible 
for the outflow. As shown by Hartmann et al. (1994), even if the initial density configuration is not toroidal, 
or magnetically-dominated, a toroidal density distribution will result from collapse. 
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strongly collimated flows even with the wide opening angles of the sheet distribution. 

We investigate the influence of ram pressure balance between wind and ambient gas by 
holding the wind speed constant and varying the ratio of the inflow mass flux to the wind 
mass flux (denoted /':) 

/' = ^. (10) 

At 

We chose to hold the wind velocity constant because young stellar object wind and jet 
velocities are well observed to be a few hundred kilometers per second (Bachiller 1996; Shu 
et al. 1991), while the range of the mass flux ratio f ' is less well known. We have focused on 
models with /' = 10,20,30,40, and 50. This was accomplished by adjusting the infall mass 
flux and fixing the outflow mass flux. We also simulated cases where the inflow mass flux 
was reduced relative to the wind, and obtained basically the same results with only minor 
differences. We note that a measure of the effect of varying /' can be seen by calculating 
the angle de at which the ram pressures are equal in our simulation. For /' = 10 we find 
de = 1.568 radians (measured from the axis) which is less than 6% of a grid cell. For /' = 50 
we find 6e = 1.519 radians which is less than 130% of a grid cell. Thus is is clear that the 
wind will be able to push all material away from the equator in the /' = 10 case while /' = 50 
simulations should experience some material attempting to cross the inner wind sphere. 



We note here the basis of our attack on this problem. The theoretical issue of the 
interaction of a stellar wind with an infalling environment is complex. There are a number 
of parameters controlling both the wind and the infalling environment. In choosing to address 
the issues involved one must ask which parameters are important enough, a-priori, to warrant 
attention? Which connect directly to issues raised in previous studies? Which allow for a 
numerically clean solution, i.e. one whose specification of initial and boundary conditions 
do not impose serious transients or artifacts on the flow? In taking on this problem we have 
adopted a strategy which, hopefully, allows to us address important aspects of the physics 
while leaving others for future works. In the present paper we focus primarily on a relatively 
simple treatment of initial/ boundary conditions (e.g. the wind and infalling environment) in 
order to clarify how a full solution to the governing equations differs from the more idealized 
solutions explored in previous works. We have constructed a focused exploration of a set of 
numerical experiments that address key issues not explored in previous studies. In future 
works we plan to open other dimension of the parameter space including a more detailed 
treatment of aspherical winds. 



3. Results 

3.1. Spherical Winds/ Aspherical Environments 

The sampling in /' described in the previous section can be thought of as a sampling 
from a relatively strong wind to a relatively weak wind, as the results will show. In addition 
to exploring different steady wind cases, the simulations serve as a crude initial exploration 
of what might be expected in models with fully time-dependent wind mass loss rates. The 
details of time- dependence are particularly important for FU Orionis stars (Hartmann & 
Kenyon 1996) where the observations show considerable variation in the wind, likely related 
to large fluctuations in the accretion rate of material through the disk. In this section, we 
start with a general overview of the results and continue by considering each simulation in 
more detail. 

The principal results can be inferred from Fig. 1 showing density snapshots of all five 
simulations. In the /' = 10 case we see the strong wind overwhelms the infall and creates a 
wide conical cavity. As /' decreases, the wind becomes weaker, the infall mass flux dominates 
and, as a result, the wind becomes more confined and collimated. Even in the /' = 50 case 
however the wind is still strong enough to avoid being completely cut off. In simulations 
where we increase the infall flux such that /' = 100, the wind is completely choked by the 
infall with its momentum slowly diffusing into the surrounding environment. Because of the 
diffusion and the direct interaction of the infall with the artificial wind boundary, we are not 
confident in the simulations after the wind becomes choked. Therefore we have only placed 
a lower limit (/' = 50) on the infall to wind flux ratio at which the wind is cut off. 

In the strong wind, (/' = 10), case (see Fig. 2), the wind pushes its way through the 
environment but is still shaped by it. By the time the swept up ambient material reaches 
the top of the grid (at ~ 180 y: see Figs. 3 and 4) it has a conical shape reminiscent of 
the outflows seen in Velusamy & Langer (1998) and Chandler et al. (1996). The outflow 
bubble opens with an angle of about 30° to the vertical. The supersonic wind and molecular 
material are each compressed at the inner and ambient shocks respectively, forming a thin 
layer about the contact discontinuity at 20° to the vertical. The dynamics of this region are 
dominated by the fact that the radially directed wind material hits the inner shock obliquely. 
Post-shock wind is redirected and forced to flow along the contact discontinuity. Eventually 
this high speed gas flows up the walls of the conical cavity and reaches a cusp at top of 
the bubble. There the focused high speed flow moves ahead of the bubble, encircling its 
nearly spherical cap with a conical "jet" . While the development of such a pattern is readily 
understood in terms of basic radiative oblique shock physics it is not clear that the cusps 
and resulting jets would be stable in real three dimensional bubbles. 
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Table 1. Model parameters. 



parameter 



value 



Simulation Domain 

Resolution 

Wind velocity, v^, 

Wind mass flux, M^ 

Central mass, M* 



Infall mass flux. Ma 



Collapse radius, tq 

Flattening parameter, 77 

Centrifugal radius. Re 



4.26 X lO"*^^ cm radius by 

8.52 X 10^^ cm along axis 

256 cells radially by 

512 cells along axis 

200 km/s 

10-^ Mq/vt 

0.21 Mq 

10^6 Me/yr x 

(/' = 10,20,30,40,50) 

5.37 X 10 

2.5 

4.28 X 10^^ cm 



16 



cm 



28.5 AU 



Fig. 1. — Comparison between bubbles at times of similar axial extent. The density of the 
/' =10,20,30,40, and 50 (left to right - strong wind to weak wind) cases are shown at their 
full heights at 100, 200, 220, 240, and 220 years, respectively. The inner wind boundary 



makes a quarter circle around the origin with radius 4.26 x 10 
aspect ratio of the bubbles as the wind becomes weaker. 



cm. 



Note the change in 
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It is interesting that as the bubble evolves, a region of warm gas (T ^ 10^ K) develops 
at the top pushing the inner and outer shocks apart there. Usually wind blown bubbles are 
classified as either adiabatic or radiative depending on whether the cooling time is longer or 
shorter than the expansion time of the bubble. We believe that this top region does not fall 
into either category, and that instead we are seeing the development of a partially radiative 
bubble. This third classification, which lies between the adiabatic and radiative cases, was 
originally explored by Koo & McKee (1992). Such a case occurs when the wind material 
resupplies thermal energy generated at the shock faster than radiative cooling can remove 
it, but when the cooling time is still shorter than the bubble expansion time. To see if the 
system is partially radiative at the top of the /' = 10 bubble we need to compare the cooling 
time tc with the wind crossing time tx- 

Given that the top of the bubble has an essentially spherical shape we can use the 
analysis of 1-D bubbles to write the cooling and crossing times in terms of the distance 
from the wind source, R, plus other simulation parameters for the /' = 10 case. A rough 
estimate of the inner shock speed is Vg = Vy, = 200 km/ s which produces a temperature of 
8.8 X IQ^K. The Dalgarno-McCray cooling curve at this temperature gives a A{n,T)/n'^ of 
1.65 X 10^^^ erg cm? / s. Using M^, = 1Q~'^ Mq/ht and the factor of 4 density increase behind 
a strong shock we can determine the cooling time from the following equation: 

t, = eth/A{n, T) ~ 4100s {R/AUf. (11) 

The crossing time, tx-, is given simply in terms of the radial distance from the origin to the 
shell and the wind velocity: 

tx = 7.5 X lO^s {R/AU). (12) 

Note that the crossing time grows linearly and the cooling time grows quadratically. Thus 
the cooling time eventually surpasses the crossing time. For the parameters given above this 
occurs when the shock is at a distance of about 182 AU or 2.7 x 10^^ cm which is well inside 
the simulation domain. We find that this is approximately the distance where the simulated 
bubble begins to exhibit its cap of warm postshock wind material. We believe that this is the 
first time that the partially radiative bubble phase has been seen in a simulation, confirming 
the prediction of Koo & McKee. 

The weaker /' = 20 case is seen in Fig. 1 to also develop a conical outer shock, but the 
overall shape is more confined or collimated than in the previous case. The cusp which forms 
at the edge of the bubble is more pronounced. This effect may due to higher inertia of the 



-11- 



Fig. 2. — The density evolution of the /' = 10 (strong wind) case. The bubble begins 
spherical but breaks out into a cone. Other features such as a partially radiative shock and 
conical "jets" appear (see text). 



shocked wind flow. At low latitudes material is channeled from regions closer to the origin 
where the density is higher. In addition, the polar regions of the bubble are expanding slower 
due to the relatively lower wind mass loss rate. Thus the refracted wind material emerges 
from the top of the bubble with higher momentum relative to the polar cap in this model 
than in the /' = 10 case, allowing the cusp to propagate farther from the cap. Another 
difference at the pole is that less warm post-shock gas exists in this case. This effect can 
also be attributed to the the slower speed of the polar sectors of the bubble. A slower outer 
shock speed implies a faster inner shock speed (in the frame of the wind). This in turn 
leads to higher post-shock temperatures and stronger cooling in the shocked wind. Note 
that the inner shock forms an angle of about 10° to vertical and the outer shock (bubble) 
opens roughly to 20°. 

At weaker winds (/' = 30) the inner shock becomes more confined and the fiow is 
even more collimated. The cusp and resulting conical "jet" is more pronounced than seen 
previously. Again this is expected because the cap of the bubble will expand at lower 
velocities. Figure 5 shows the details of the fiow pattern for this case. Since the /' = 20 
and /' = 30 models are similar the figure emphasizes the points concerning the fiow pattern 
made above. The wind shock opens about 2° to vertical and the bubble opens about 20°. 

In the second weakest case (/' = 40) the inner shock becomes almost elliptical or 
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Ambient Shock 



Contact Discontinuity (CD) 




Cusps ^ 



Fig. 3. — A cartoon of the /' = 10 (strong wind) case, labeling the location of the important 
fluid-dynamics features mentioned in the text. The central source and freely expanding wind 
are shown. The major surfaces of discontinuity (shock, contact) are identified as well as the 
sharp cusps in the wind shock. The flow of the redirected wind is shown schematically. 
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Fig. 4. — The density, pressure, temperature, and velocity of the /' = 10 (strong wind) case 
at 180 years. The velocity magnitude is represented by grayscale shading and the direction 
is represented by overlaid normalized arrows. Note that in the freely expanding wind region 
(white) the velocity vectors have been sampled at regular grid points and fool the eye into 
seeing a precollimated wind. The wind in that region is actually expanding spherically. 
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Fig. 5. — The density, pressure, temperature, and velocity of the /' = 30 (moderate wind) 
case at 220 years. 
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Fig. 6. — The density evolution of the /' = 50 (weak wind) case. Confinement is strong, and 
is maintained by a prolate shock. 



"bullet" shaped. The tip of the shock at the axis is suspect because of the imposition of 
refiective boundary conditions and may not form if the cylindrical symmetry is relaxed (see 
for example Stone & Norman 1994). The wind is strongly focused by the oblique inner shock 
producing an almost vertically directed high velocity fiow. The resulting bubble is very well 
collimated. In this case the conical cusp is not as pronounced and the dynamics look very 
much like a jet. At the tip of the outfiow (along the axis) the the inner shock is acting in the 
same manner as a jet shock, decelerating vertically directed wind and shunting it off in the 
direction transverse to the propagation of the outfiow. The full opening angle of the outer 
shock of the bubble is about 35° (17.5° to vertical). 

The weakest wind case (/' = 50) is shown in Figs. 6, 7. It exhibits very strong 
collimation. The inner shock is strongly prolate and closes back on itself at a relatively small 
distance from the central source, and redirects all the wind into the collimated outfiow. There 
is no extended region of freely expanding wind. This occurs because infiowing material is 
overwhelming the wind at low latitudes, inhibiting its expansion everywhere except at the 
poles. Even though the opening angle of the bubble's outer shock at the base is 35°, the 
resulting outfiow will appear cylindrical because of the focusing. 

The trend of the evolutionary timescales of the different simulations allows us to infer 
the dynamical consequences of collimation. The timescales are taken from the snapshots 
in Fig. 1 which were chosen such that the bubbles are of comparable height. As the wind 
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Fig. 7. — The density, pressure, temperature, and velocity of the /' = 50 (weak wind) case 
at 220 years. Note the focusing effect of the inner shock evident in the velocity plot, and the 
wide region of vertically directed flow. 
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becomes weaker (/' = 10 to /' = 50) the timescales go 100 y, 200 y, 220 y, 240 y, 220 y. The 
bubbles decelerate as the wind weakens. The effect saturates however and then reverses as 
the colhmation becomes stronger. At first glance this reversal may be unexpected because 
as winds become weaker they should supply less momentum to the bubble. However, as the 
wind becomes focused into a jet the net flux of momentum, F = pv'^A, is injected into the 
environment across a smaller angular extent A. Thus the tip of a highly collimated bubble 
expands at a higher velocity than would occur for a spherical bubble. 

We shifted from the "strong wind" to "weak wind" cases in these simulations by increas- 
ing the infall mass flux while flxing the wind velocity. If instead we shift from /' = 10 to the 
/' = 50 by decreasing the mass flux of the wind and keeping the infall mass flux flxed we 
flnd that the results are qualitatively the same. The post inner shock region expands slightly 
faster due to the lower density and reduced cooling there, but the effect is not signiflcant. 
The results are, therefore, determined solely by the parameter /'. 

The results of the simulations show the complicated way in which the environment 
shapes the outflow. This occurs directly at the outer shock by providing an inertial gradient 
along the shock face. It also occurs indirectly by affecting the shape of the inner shock. 
Shaping the inner shock leads to shock focusing of wind material, creates cusps at the tips 
of the bubble and, in the most extreme cases, leads to focusing into a jet through prolate 
inner shocks. In all cases, the simulations show the "internal" flow of post inner shock wind 
material plays an important role in the dynamics. 

The results also show the way in which outflows affect the environment. The winds 
carve out evacuated cavities of potentially large angular extent and sweep ambient material 
at a high rate. The swept up mass totals 1 x 10~^MQ/yr to a few times 10~'^Mi7)/yr as the 
wind becomes weaker. Because the winds become more collimated and on the whole expand 
much more slowly as they weaken, ambient mass is swept at a slower rate. The wider angle 
winds are able to collect mass much more rapidly even though the envelope has less material 
in it. 



3.2. Spherical Winds/Aspherical Environments 

Although they are not the focus of this study, we have included the results of simulations 
with aspherical winds injected into a spherical environment for comparison. We note that 
addressing the issue of aspherical winds means specifying a model for those winds. This 
raised a number of questions suh as; should simple sinusoidal variation be used or something 
akin to the cylindrical stratification seen in wide-angle models like those of the X-winds? If 
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we focus on the former should we vary only the density from pole to equator or only the 
velocity. If they both vary then should momentum or energy flux be held constant across 
the stellar surface? Given these the complexities we have chosen a highly simplified set of 
boundary conditions for these models and present them as means of contrasting the results 
of the previous section. 

The distribution in the environment is given by Terebey et al. (1984) with the same 
infall parameters as used in the strong-wind simulation described above: essentially it is the 
previously described sheet distribution except unflattened. The wind velocity distribution is 
given by 

v^(cos^)^°s("°^^=)/'°35, (13) 

where 9 is the polar angle, v^ is the maximum polar 200 km/s velocity and the angle of 
half-maximum velocity is 9c. The wind mass density at every point is kept the same as in 
the strong-wind simulation. 

The density of the result of injecting a mildly focused wind {9c = 70°) is shown in 
Figure 8. Even at 300 years, the bubble is much more limited in size than in the sheet 
environment cases: the expansion speed is lower and very little warm gas exists to push the 
bubble outwards. The focusing is also much more limited because of the the environment. 

The resulting density snapshot of injecting a more focused wind {9c = 30°) into the 
spherical environment is shown in Figure 9. Here the focusing is pronounced, but the 
bubble extent at 300 years is still comparatively small because of the greater mass in the 
polar regions. As is to be expected, a focused wind in a spherical environment can also 
produce a focused outflow. Detailed results, however, such as the precise bubble shape 
should not be taken as final because the infall material is able to press up against the 
artificial wind boundary. 

These simulations show that the shape of a wind blown bubble is not unique to a single 
set of boundary conditions. As one might expect we do see coUimated outfiows emerging from 
collimated winds though there are some differences from the spherical wind case in terms of 
propagation speed and post-shock temperatures. These results show that the issue is always 
one of wind momentum or ram pressure injection as a function of angle. With spherical- 
winds/aspherical-infall shock focusing helps redirect the wind ram pressure towards the pole. 
This is explored in more detail in the next section. In the aspherical winds/spherical-infall 
case the wind momentum must be strongly collimated a-priori to produce a narrow outflow. 

Thus flow shapes are not unique nor should one expect them to be. However, as we 
will explore in the next section, the simulations presented in section 3.1 demonstrate an 
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Fig. 8. — The density at 300 years of a case where 
a non-spherical wind is injected into an isotropic 
infalhng environment. The half- maximum veloc- 
ity is at 6c = 70° from the pole. Note the small 
extent of the bubble at this late time and the mi- 
nor focusing. 



effectiveness of wide-angle wind models not seen in previous studies. 



4. Discussion 

In this section we discuss the results of our simulations in light of previous analytical 
models of YSO wind blown bubbles We also attempt to establish some points of contact 
between our models and molecular outflow observations. In making these comparisons we 
must note a limitation of the present simulations, namely that they formally refer to only 
very small spatial and temporal scales. The overwhelming majority of observed molecular 
outflow structures have typical lifetimes of 10^ — 10^ yr, in contrast to the 10^ yr time 
sequences shown here. In addition, for flows that are much larger than 0.1 pc, it is very 
likely that the bulk of the swept-up interstellar material originally resided in regions outside 
the protostellar core, let alone the inner infall regions. Nevertheless, there are a couple of 
observations of younger and smaller outflows (Chandler et al. 1996) where our simulations 
may be more directly relevant. In addition some of the physical results may be generalized 
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Fig. 9. — The density at 300 years of a case 
where a narrow wind is injected into an isotropic 
infalhng environment. The half-niaxiniuni veloc- 
ity is at Oc = 30° from the pole. The outflow 
is coUimated as expected. Here, as in the wider 
wind case in Figure 8, we find that the bubble 
is growing much more slowly than those formed 
in the collapsing sheet environment. The inf ailing 
material is reaching the wind sphere in this case, 
so the exact shape of the bubble obtained in this 
simulation is suspect. 



to larger-scale situations. 



4.1. Comments on the Shu, et. al. "Snowplow" Outflow Model 

In Shu et al. (1991) a simple, elegant "snowplow" model of the formation of bipolar 
molecular outflows is presented. The basic features of this model were a central wind with 
a mass flux varying with polar angle and an environment with a density that varies with 
polar angle and falls off in radius as 1/r^. A number of simplifying assumptions are invoked. 
Cooling is assumed instantaneous so that the bubble is a purely momentum-driven thin shell. 
The interaction of the ambient material and the wind is taken to be fully mixed and ballistic 
such that the dynamics at a given polar angle are independent of those at nearby angles. In 
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addition, the swept up wind mass is taken to be negligible. 

In the snowplow model, changing the ratio of total infall to wind mass flux (/') only 
effects the bubble's expansion timescale but doesn't change the self-similar shape. The 1/r^ 
ambient profile leads to a simple expression for the bubble shape and the variation with polar 
angle in the extent of the bubble is purely determined by the angular variations of wind and 
ambient medium. The simulations however show a strong variation of bubble shape with /'. 
This discrepancy implies that some dynamical aspects of the simulations are not captured 
in the snowplow model. 

The most important physical difference between the simulations and the snowplow model 
is that the snowplow model does not account for the post inner shock redirection of wind 
material. Pressure gradients in the post-shock shell that come about because of the variation 
in the ambient density, and resulting variations in shock obliquity and shock speed have no 
consequence in the snowplow model, but have profound effect in the simulations. If we 
were to somehow model the non-radial direction of the postshock wind and use that in the 
snowplow model we might obtain a better agreement. The question becomes should one use 
the flow inside or outside the inner wind shock as the input wind to the snowplow model. 
While this may seem, on one level, to be a question of where one starts the calculation, 
more fundamentally it is an issue of what initial physics one includes to determine the wind 
distribution. Li & Shu (1996) make a modiflcation of the snowplow model in this way by 
incorporating the distribution of the wind after it has been accelerated and partially focused 
by magnetic launching. 

If dynamically important physics in wide wind scenarios are left out of the snowplow 
model, it leads to the question whether the model is useful in predicting observations. Mas- 
son & Chernin (1992) and Li & Shu (1996) have used the snowplow model to calculate 
distributions of mass and velocity to argue whether or not wide angle winds can drive out- 
flows. Our simulations, however, show that even if the driver is initially a wide angle wind 
it can be focused by a nonlinear interaction between the wind and environment. Our simu- 
lations appear to blur the distinction between wide angle winds and jets as the sole driving 
mechanism of the outflow. 



4.2. Mixing and momentum transfer 

A strong assumption used in the analytic models of Shu et al. (1991) and Wilkin & 
Stahler (1998) is that postshock wind and postshock ambient material become fully mixed 
on a dynamically short timescale. This allows the shell to be treated as a single fluid and 
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greatly simplifies the calculations. To support the assumption, Wilkin & Stahler consider 
supersonic shear flows within the postshock regions and state that some form of instability 
will generate turbulence and rapidly mix the wind and ambient fluids. Our simulations have 
some diffusive mixing of fluids due to numerical effects, but no complete mixing occurs and 
there is no turbulence modeled in the code. Since the simulations show the post-shock shear 
flow influencing the shape of the bubble, it is worth discussing when turbulent mixing of 
ambient and wind material will occur and, if it does occur, how it might affect the bubble's 
evolution. 



4.2.1. Turbulent length scale 

The onset of turbulence in jets and bubbles is still poorly understood, but one estimate 
of the momentum transfer rate between material in a jet and a stationary ambient medium 
has been performed by Canto & Raga (1991). From this, they calculate the length, Lf, a 
laminar jet can travel before it is subsumed in a turbulent boundary layer flnding that Lt 
scales linearly with jet Mach number, and inversely with the jet radius Rj. At Mach 10 
they flnd Lt/ Rj ~ 170. The calculations were performed assuming slab symmetry and are 
applicable to our supersonic postshock wind regions. Raga et al. (1995) extend the results to 
the case where both fluids are moving and flnd that the spreading rate is actually reduced. 

Using the result of Canto and Raga we can estimate when the flows inside the bubble 
become turbulent. For this purpose we calculate the width of the shell of post-shock wind in a 
slightly aspherical isothermal wind-blown bubble. The calculation is simplifled by the nearly 
spherical shape, but is still consistent with shock focusing because even mildly aspherical 
bubbles will produce focused flows (Frank & Mellema 1996). If density p^^ in the postshock 
region is constant, the mass in that region, M^g, can be written in terms of the volume V^s 
as: 

M^s = V^s Pws- (14) 

If the speed of the shock Vg is small compared with that of the wind, M^s can also be written 

M„, = M^ t, (15) 

where M^ is the mass flux of the wind and t is the age of the bubble. If the shock is 
relatively thin and the bubble is not too aspherical then we can approximate the volume 
as Kos ~ AirRgAR^s, where Rg is the shock radius and AR^g is the shock width. In terms 
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of the shock Mach number M and the preshock wind density pw the isothermal condition 
requires that p^^ = pwM'^ = M^M^/i'inR'^^v^) (Shu 1992). Equating (14) and (15) and 
solving for the shock width gives: 

Afi™ = l^. (16) 

Once we calculate a shock speed, we can calculate the relative width of the bubble, ^R^sl Rs- 
For simplicity we choose the ambient environment to fall as 1/r^, which, as mentioned in the 



previous section, implies a constant shock velocity of v^ = \J [M^^/ Ma)vu, (a/2). Choosing 
typical values of Vw = 100 km/s, M^/Ma = 1/10, M = 100, a = 0.2 km/s, we obtain: 



AK„,, / M,„ a 



Rs {Ma'i M2 100 



;i7) 



Note the high Mach number is a result of the fact that in an aspherical bubble the shock 
creating the shear flow will be oblique and therefore will not strongly decelerate the wind 
material. In addition the Mach number is an inverse function of the sound speed. The 
strong cooling in the postshock region keeps this speed low, which also helps to keep the 
Mach number high. The simulated shear flow has velocity on order Mach 10 (as opposed to 
the compression Mach number M = 100) and we can couple equation (17) with the result 
of Canto and Raga, to flnd that 

^^^r^l. (18) 

Rs AR ^ ' 

This suggests that the flow becomes turbulent only on a size comparable to that of the 
bubble, and that the focused tangential flow will have time to influence the dynamics before 
mixing occurs. 

Two features this calculation doesn't take into account further reduce the chance that 
turbulence will dominate the dynamics. First, the continuous supply of fresh wind material 
along the length of the shear flow is ignored as it was in the Canto & Raga estimate. This 
effect would allow unmixed wind material to be continually resupplied at high latitudes, 
possibly at a rate faster than the turbulence incorporates material. Second, magnetic flelds 
are ignored. Fields may thread the shell especially when the bubble is small. There is some 
evidence that toroidal flelds, which are present in magnetocentrifugal wind launching models, 
retard the entrainment of envelope material (Rosen et al. 1999). Poloidal flelds at moderate 
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Mach numbers can also become wrapped up in vortices in the mixing fluids, strengthening 
the field which in turn stabilizes the shear fiow (Frank et al. 1996). 

We have neglected effects in our calculation that may play a mediating role. The most 
important lie in the crude way in which we treat the geometry of the non-spherical shock. We 
also haven't considered the early transition from a spherical to an elliptical bubble, when the 
shear fiow Mach number will be more modest. Also, when the shear is stronger, centrifugal 
forces in the bubble may be important. The details of what instabilities form and how they 
grow in these geometries could affect the growth rate of the turbulence. Thus while the 
estimate (18) is suggestive, how quickly the flow becomes turbulent merits further study. 



4-2.2. Direction of the turbulent flow 

Taking the opposite view, we can assume that turbulence will be important dynamically 
at some time. How then would the resulting mixed fluid behave? Wilkin & Stabler (1998) 
assume the postshock flow is dominated by the shocked ambient gas and will be carried down 
towards the disk. While there is much more mass in the ambient material, the momentum 
in the wind material is quite high and the flow could be driven poleward. Below we calculate 
the direction a mixed flow would take in a wind blown bubble shell using a simple model. 

In flgure 10 we show preshock and postshock regions for the ambient and wind material 
with a mixing layer sandwiched between. These regions denoted by subscripts a,w, and L 
respectively. We assume the infall Va and wind velocities Vu, are constant and the turbulent 
boundary layer grows into the postshock regions at the same constant velocity Vq. In this 
framework the question becomes: is vl positive (i.e. poleward) or negative (i.e.equatorward)? 

Mass conservation implies: 

Pl2Vo = psaVo - pswi-Vo) ^ pL = -{psw + Psa)- (19) 

Likewise, momentum conservation for the component moving parallel to the faces of the 
boundaries implies 



{PlVl)2Vo = {psaVsa)Vo - ipswVsw)i-Vo) => ^PlVl = psaVsa + PswVsw (20) 

The preshock densities are determined from the mass flux and velocity through p = M/(47ri?^f ) 
assuming the shell is thin. By taking the shocks to be isothermal, we obtain the postshock 
densities by multiplying the preshock densities by M^ where M is the Mach number. Because 



-25- 



Wind 




Pw = 



M 



w 



4^R|v^ 



Sliocked 
: Wind 
IVIaterial 






L 


■.:■: 




-.■.■.- 






^ 


. ■ ■ ' 


-vo 


■'■'.■■'■ 




■ ■. '■ 


Ps 


~ 2 


■■:'-'. 



Turbulent ;J 
Boundary ;g 
Layer 5>; 



':-%«« 



■vlIII 



Pl? 



Shocked 
Ambient 
Material 



V^ 



-VgSinOsa 



Psa = PwMsa 



Ambient 
Medium 



V„ 



'sa 



Pa = 



M, 



47tR|Vg 



Fig. 10. — Cartoon for a mixing model. The diagram represents a small section of the 
bubble's shell with wind and ambient material encountering shock layers at an oblique angle. 
The central mixing layer subsumes material from postshock wind and ambient material which 
both have velocities parallel to the layer interfaces. The text describes the calculation of the 
direction of the velocity in the mixing layer. For parameters of interest, we find that the 
wind determines the direction of the mean velocity of the turbulent region. 
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the fluid velocity parallel to the shock remains unchanged through the front, Vsa = —Va sin 6sa, 
and Vsw = VyuSinOgw Incorporating these velocities and solving the equations (19) and (20) 
for vl gives: 

We are most interested in whether the mixed flow moves in the direction of the wind or 
the infall. If the wind and ambient velocities are mostly aligned across the shell, then 
sinOsa = sin ^5^, and the condition for the mixed flow moving in the direction of the wind is: 



M^ / M,^ 



Ma \Msa 



> 1. (22) 



For typical values of Ms^, ~ 100, Msa ~ 10, M^/Ma = 1/10, the left hand side of (22) is 
~ 10. Thus, the flow moves in the direction of the wind. 

We find that the large momentum input of the wind would tend to accelerate the mixed 
material poleward. At the pole, material collides and sprays ahead of the bubble. This is the 
Canto focusing mechanism Canto & Rodriguez (1980), called conical converging flows, which 
operates in simulations under a variety of circumstances (Mellema & Frank 1997). Material 
would be cleared and accelerated above the pole leading to a more elongated bubble, and 
possibly a very narrow jet. 

Ignoring shear flows simplifies calculations greatly. We have argued, however, they are 
important dynamically in bubbles driven into stratified environments whether turbulence is 
dominant or not. Future models will have to consider such flows in more detail. 



4.3. Properties of young outflows 

As noted in the introduction. Chandler et al. (1996) point to three qualities that models 
of molecular outflows must reproduce to match their observations of TMC-1 and TMC-lA. 
They require conical outflow lobes close to the star, evacuated outflow cavities, and moderate 
30° — 40° opening angles, and they argue that existing models do not explain these effects. 

We flnd, however, that our simulations create outflows meeting the requirements. The 
opening angle of bubbles in the /' = 20 to /' = 50 cases are 35° — 40°. To be sure, these 
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angles are determined from the density profiles of the simulations and should in the end be 
determined from synthetic emission maps. However, the angles are not as wide as might 
naively be expected from a wide angle wind. The similarity to TMC-1, TMC-IA and other 
young outflows is encouraging because these simulations are on a small spatial and time scale 
(300 years, 10^^ cm) and would hopefully match up more readily with younger outflows such 
as these. 

The opening angle of the bubble in the /' = 10 case is 60°, which is too wide to compare 
with TMC-1 and TMC-lA, but wider outflows have been observed. The outflow around IRSl 
in Barnard 5 has an opening angle of about 125° (Velusamy & Langer 1998). Although these 
outflows may be explained by multiple non-aligned driving winds or jets, very strong wide 
angle winds may be at work in such situations. 



4.4. Mass vs. velocity relations. 

Masson & Chernin (1992) used CO line intensities of NGC 2071, L1551, and HH 46-47 
and assumptions about optical depth to obtain the amount of mass per velocity channel 
in these outflows. They measured a power law in this mass versus velocity relationship 
of Am/Av = v'^ with an exponent of 7 = —1.8. With such a steep slope little mass 
is accelerated to high velocities, leading Chernin and Masson to conclude that any high 
velocity driver of the molecular outflows must have a small cross-section and consequently 
cannot be a wide angle wind. Calculations of the observed mass vs. velocity relation in 
other outflows have been calculated with similar results but with interesting differences. 
Chandler et al. (1996) have pointed out that it is hard to justify fitting the distributions in 
their observations with a single power law, although they do try to fit them with two. They 
find that some "average" power law curve obtains a 7 ~ —1.8. Bally et al. (1999) bring up 
the issue of sensitivity of the transformation from line intensities to masses to assumptions 
of optical depth. The universality of a mass vs. velocity relationship may therefore be in 
question, whether the relationship is even a power law at all, and whether the mass has 
been properly calculated. We are interested however in making some comparison of models 
and observations, so it is worth discussing the mass-velocity relation calculated from the 
simulations. 

Figure 11 shows the histograms of ambient mass as a function of projected velocity from 
our simulations. We able to unambiguously identify the ambient material through the code's 
fluid tracking and independently by the location of the contact discontinuity. The swept up 
masses are 1-2 orders of magnitude lower than the masses shown in the TMCl and TMCIA 
observations of Chandler et al. (1996). However, our numerical simulations are limited in the 
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Fig. 11. — Histograms of mass per line-of-sight (LOS) velocity channel for the simulations 
at different times. As the simulations evolve they sweep up more mass, thus the older curves 
are above the younger curves. The velocity channels are 0.6 km/s wide and the LOS is 45°. 
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length of time considered and the spatial region that can be encompassed while maintaining 
adequate resolution. Thus, our simulations correspond to sweeping up much less envelope 
material than has occurred in most observed outflow regions. Quantitatively, our models 
correspond to extremely young, rare systems; qualitatively, many of the effects we show 
should be present in older objects. 

As found in the Chandler et al. (1996) observations the simulation curves are not well fit 
by a single power law. The curves are logically split into two parts: a "shallow" portion that 
can be fit on the simulation plots roughly by a power law of exponent ranging from —1.5 to 
— 1.3 and a steeper portion. The shallow portion matches well with the TMCl and TMCIA 
results of -1.26 and -0.75 respectively, and is fairly close to the Chernin & Masson value of 
7 ~ —1.8. These fits are subject to some small variation from the choice of where one portion 
begins and ends, and the size of the velocity bins, but varying among reasonable choices for 
these only changes the exponents by a tenth or so. The steeper portions of the simulation 
curves are fit with exponents steeper than —10, to be compared with an observation of ~ —3 
(TMCl and TMClA) and a jet-model value of —5.6 (Zhang & Zheng 1997). There appears, 
however, to be some fiattening of the steep portion as the system evolves in the /' = 50 case, 
which is the most focused model. The slope fiattens from —19 to about —7 to —3.25 by 220 
years. The focusing accelerates the polar tip and sweeps up increasingly more mass to higher 
velocities. This acceleration might also occur in other less focused cases when the outfiow 
"breaks out" of the protostar's cloud core into much less dense material. If this breakout 
fiattens the mass vs. velocity curve at high velocity then it could explain the observationally 
derived curves from older outfiows. 

The shallow and steep portions of the mass-velocity curves are connected by a "hump," 
which does not appear in the observationally determined curves. This feature is also seen 
in the mass-velocity plot of Zhang & Zheng (1997) and we find that in our simulations it is 
largely due to the innermost portions of the ambient shock - those farthest from an observer. 
Optical depth effects could play a role in translating this part of the curve to an accurate 
integrated CO line profile. 

The mass vs. velocity relations for the simulations look promising, with one caveat 
being that questions remain about this method of matching observations to models. Another 
caveat is that extrapolating these results to those of large-scale outfiows (larger than 0.1 pc) 
is not straightforward. Large outfiows are likely to involve significantly more material from 
regions exterior to the cloud than material in the infall in the infall region, and thus the 
simulations here may refer only to structures on 100 - 1000 AU scales. The Chandler et. al. 
observations are more relevant in this case. Noting these provisos, the models have mass vs. 
velocity relations whose shallow portions are about as steep as those quoted by Chernin & 
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Masson, and the steeper portions may flatten as the flow focuses or breaks out of the cloud 
core. The simulated relations have striking similarities to the observationally derived curves 
despite the fact that the model systems are driven by wide angle winds. 



5. Conclusion 

We have performed radiative simulations of the interaction of a central spherically sym- 
metric wind with an infalling environment under conditions similar to those found around 
young stellar objects. We include gravity from the protostar and rotation about the axis. We 
have examined the effect of varying the relative infall to wind mass flux on the morphology 
and kinematics of the resulting outflow. 

In our simulations, shaping by the toroidal YSO environment of the wind shock, in 
addition to the outer molecular shock, focuses the outflow, leading to morphologies and 
kinematics not explored in earlier analytical models. Most importantly, the outflows become 
progressively more collimated as the central wind becomes weaker relative to the infall. 
Stronger winds can reverse the inflow of material creating the wide evacuated conical cavities 
seen in some observations of molecular outflows. Weaker winds are focused more tightly and 
form elongated, bipolar structures. The simplifled snowplow model for these outflows does 
not explain these variations. Thus, conclusions made from the snowplow model that wide 
angle winds cannot drive molecular outflows should be reexamined. 

We suggest that the effects of turbulent mixing of wind and ambient material, which are 
not treated explicitly in the simulations, will not signiflcantly alter the results. The mixing 
is shown to be a relatively inefficient process. The shear speed of the post-shock wind in the 
simulations is high and we estimate that the poleward focused flow will traverse a large part 
of the wind blown bubble's circumference before mixing becomes signiflcant. If the mixing 
does occur, the momentum in the wind is large enough to force the flow towards the poles 
and possibly form a conical converging flow. 

We find a favorable comparison of the mass vs. velocity relations derived from the 
simulations to those from observations. The shallow portion of the curve at low velocity can 
be roughly fit by a power law of exponent from about —1.5 to —1.3. The steeper portion 
at high velocity is fit by an exponent of —10 or steeper, although the most focused of our 
simulations shows a gradual fiattening. This leads us to speculate that when the other 
simulations run to breakout of the cloud core there might be further focusing and a similar 
transition in the mass vs. velocity curve. 

These simulations show that directed outfiows can result from wide angle winds. The 
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nonlinear interaction of wind and environment leads to dynamics that are not easily modeled 
outside of simulation, but which lead to behavior similar to what is seen in young outflows. 
Properties that have been previously been proscribed to jet driven systems can be found in 
the wide angle winds. Further synthetic observations need to be performed on simulations 
like those in this paper to indicate the other ways in which the wide angle wind model offers 
an explanation for molecular outflows. 

For low mass YSOs it is probable that wide angle winds will be generated by some form 
of MHD mechanism associated with an accretion disk (Ouyed & Pudritz 1997) or disk-star 
boundary layer (Shu et al. 1994). Although our simulations are hydrodynamic and involve a 
maximally uncollimated (i.e. isotropic) wind, the shaping of outflows should be general and 
occur even in the presence of magnetic fields. This is because shock focusing could occur 
even when the shocks are magnetized. In addition, MHD driven wide angle wind models 
will likely enhance the outflow collimation in two ways. The first is because these models 
typically collimate the wind early, in the "wind sphere" region we do not model, to send 
more momentum along the poles. The second is because toroidal fields are generated and 
carried out by the wind in such models, providing radially directed hoop stresses in the wind 
post-shock region. 

In future papers we hope to explore in more detail the emissive properties of wide wind 
driven simulated outflows and make a more direct comparison with observations, as well as 
study long term behavior. Additionally, these studies of the relative magnitude of the wind 
have laid the foundation for an examination into how time variation in the wind affects the 
outflow. With time dependence of the source, which we expect on observational grounds, the 
outflows can be extremely complicated kinematically. Finally a more detailed exploration 
of the parameters associated with the asphericity of the wind should be explored allowing a 
better link between its shape and the resulting properties of the outflows. 

We wish to thank Thomas Gardiner for his helpful discussions. This work was sup- 
ported by NSF Grant AST-0978765 and the University of Rochester's Laboratory for Laser 
Energetics through a Frank J. Horton Fellowship and their computing facilities. 



A. Collapsing Sheet Model 

In this appendix we briefly recapitulate the equations of the sheet density distribution 
described in Hartmann et al. (1996). 

The collapsing sheet is described as a function of spherical radius from the origin r and 
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cosine of the polar angle /x. In spherical coordinates the mass density and velocity is: 



-(-,/i) = (sech(r^^o)r(-^)^^^(^,/i) (Al) 

Pn Re \tanh(?7)y pn Re 

pn Re \ReJ \ Po J \P0 Jf^ J 

-(f,/^) = -(^r'^"(^ + -Y' (A3) 

Vk Re Re \ PoJ 

^(f,.) = (^)-"H,„.,)(^l±£-y (A4) 

Vk Re Re Va*o(i -/i^)/ 

Wfc Re Re V A*o/ VI -/^ / 

where /xq = Po{r / Re, p) specifies the location (cosine of polar angle) on a reference sphere at 
some radius Tq which a gas parcel at the given coordinates came from. It is defined implicitly 
by: 



r 1-fil 



Re 1 



(A6) 



The parameter r] is the ratio of the distance inside-out collapse has progressed into the 
sheet To, to the scale height H of the initial, static, self-gravitating sheet. It also describes 
the degree of fiattening in the central density distribution. The value of rj is taken as 
constant during the simulations. Other basic constants include the centrifugal radius. Re, 
the Keplerian velocity at the centrifugal radius, Vk, and a density p„. These in turn can 
be defined using the mass of the (forming) central star M*, and an angular velocity at the 
radius r^ of 0o- 



R^ = ^ (A7) 



p = M^L (A9) 
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